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given  the  task  of  detecting  the  signals  on  data  from  a  single 
instrument.  The  false  alarm  rate  was  set  at  3  to  4  per  hour  by 
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detector,  11.  This  relative  ranking  is  in  agreement  with  an  earlier 
study  using  signals  buried  in  phase-scrambled  noise. 

The  standard  detector  detected  14  out  of  the  15,  whereas  it  was 
worse  than  all  the  other  three  in  the  early  study.  This  may  be  because 
the  standard  detector  takes  the  absolute  value  of  the  filtered  data 
instead  of  squaring  it  and  thus  is  more  resistant  to  the  effects  of 
spikes  and  weak  precursor  signals.  More  likely  the  standard  detector 
was  at  a  disadvantage  in  the  earlier  study  becuase  the  signal  time 
window  was  held  at  2,5  seconds  while  the  signal  windows  of  the  other 
detectors  were  much  larger  as  appropriate  to  the  larger  signals  which 
were  buried  in  noise. 

All  existing  detectors' are  probably  inadequate,  in  their  present 
form,  to  be  reliably  used  to  select  data  segments  for  transmission 
because  of  weaknesses  in  the  rejection  of  spikes,  drop-outs,  and  noise 
bursts;  weakness  in  detection  in  coda;  and  inability  to  suppress 
multiple  detections  in  a  strong  event.  However,  all  detectors  are 
equally  susceptible  to  improvement  in  these  aspects. 
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ABSTRACT 


Three  advanced  detectors  and  a  standard  detector  were  tested  on  15 
of  the  weakest  40  teleseismic  signals  associated  to  events  and  received 
at  an  array  in  Alaska  in  the  time  interval  1-6  November  1980.  The 
signals  were  originally  detected  on  array  beams  and  the  detectors  were 
given  the  task  of  detecting  the  signals  on  data  from  a  single 
instrument.  The  false  alarm  rate  was  set  at  3  to  4  per  hour  by 
comparing  single  instrument  detections  to  beam  detections.  The  optimum 
2-1 og  detector  detected  all  15;  the  Walsh  detector,  13;  and  the  MARS 
detector,  11.  This  relative  ranking  is  in  agreement  with  an  earlier 
study  using  signals  buried  in  phase-scrambled  noise. 

The  standard  detector  detected  14  out  of  the  15,  whereas  it  was 
worse  than  all  the  other  three  in  the  early  study.  This  may  be  because 
the  standard  detector  takes  the  absolute  value  of  the  filtered  data 
instead  of  squaring  it  and  thus  is  more  resistant  to  the  effects  of 
spikes  and  weak  precursor  signals.  More  likely  the  standard  detector 
was  at  a  disadvantage  in  the  earlier  study  becuase  the  signal  time 
window  was  held  at  2.5  seconds  while  the  signal  windows  of  the  other 
detectors  were  much  larger  as  appropriate  to  the  larger  signals  which 
were  buried  in  noise. 

All  existing  detectors  are  probably  inadequate,  in  their  present 
form,  to  be  reliably  used  to  select  data  segments  for  transmission 
because  of  weaknesses  in  the  rejection  of  spikes,  drop-outs,  and  noise 
bursts;  weakness  in  detection  in  coda;  and  inability  to  suppress 
multiple  detections  in  a  strong  event.  However,  all  detectors  are 
equally  susceptible  to  improvement  in  these  aspects. 


-2- 


TABLE  OF  CONTENTS 


Page 

ABSTRACT  2 

LIST  OF  FIGURES  4 

INTRODUCTION  5 

EXPERIMENTAL  DESIGN  10 

RESULTS  AND  DISCUSSION  12 

REFERENCES  14 

APPENDICES  1-1 

I.  April  7,  1981  Geotech  Memorandum  from  R.  R.  Blandford,  1-1 
"Evaluation  of  SDAC  On-Line  Detector  on  Nov.  1,  2, 

1980  -  On-Line  ALK  Data  _  Revision  of  30  March  1981 
Memo 

II.  Program  Descriptions  and  Condensed  Listings  II-l 

Z  Detector  II-2 

Walsh  Detector  II-9 

MARS  Detector  11-14 

On-Line  Detector  11-18 


I 

I 


LIST  OF  FIGURES 


Title 

False  alarms  versus  relative  magnitude 
threshold  for  different  detectors  on  Pinedale 
test  tape. 

False  alarms  versus  relative  magnitude 
threshold  for  different  detectors  on  NORSAR 
test  tape. 


INTRODUCTION 


In  a  recent  study,  Blandford,  Racine,  and  Romine  (1981)  showed  that 
the  best  detector  of  signals  buried  in  phase-scrambled  noise  was  one  in 
which  first  an  optimum  S/N  filter  was  updated  and  applied  to  the  data. 
The  filter  was  followed  by  computation  of  a  10-second  average  of  the 
squared  filtered  data;  the  logarithm  of  this  short-term-average  (STA) 
was  taken  and  the  result  used  to  recursively  update  a  long-term-aver age 
(LTA)  with  a  time  constant  of  2  minutes.  The  standard  deviation  (o)  of 
the  difference  between  the  STA  and  LTA  was  also  recursively  computed  and 
finally,  a  Z-detection  statistic  (STA  -  LTA) / a  computed.  On  Figures  1 
and  2,  which  are  taken  from  the  study  by  Blandford  et  al . ,  we  see  that 
this  detector,  indicated  by  the  letter  Z,  does  have  the  highest  relative 
magnitude  threshold  for  a  fixed  false  alarm  rate. 

The  next  best  detector,  on  the  average,  was  the  Walsh  detector, 
designated  in  these  figures  by  the  letter  W.  This  detector  is  described 
in  detail  by  Goforth  and  Herrin  (1981).  The  data  is  prefiltered  with  a 
recursive  frequency  filter;  the  Walsh  transform  is  taken  of  successive 
3.2  second  windows  with  1.6  second  overlap.  The  absolute  values  of  the 
Walsh  coefficients  corresponding  to  the  range  of  circular  frequencies 
with  good  S/N  are  summed,  after  each  is  weighted  with  its  own  long-term 
average.  This  sum  is  the  detection  statistic,  and  the  threshold  is  set 
by  maintaining  a  histogram  of  512  such  successive  values  (a  time 
interval  of  about  13  minutes)  and  by  comparing  each  value  to  a  constant 
times  the  difference  between  the  50th  and  75th  percentile  histogram 
entries.  If  the  threshold  is  crossed  twice  consecutively,  a  detection 
is  declared. 

Less  than  0.01  m,  units  worse  than  the  Walsh  detector  are  the 

D 

"standard"  detectors  (run  18  for  Pinedale  and  24  for  NORSAR) ,  a  simple 
upgrade  of  the  online  detector  in  which  the  short  term  average  is 
squared  instead  of  the  absolute  value  being  taken,  and  in  which  a  single 
threshold  crossing  is  sufficient  for  detection  instead  of  2/3  or  3/3. 
These  changes  are  in  accord  with  theories  of  detection  as  discussed  by 
Blandford  et  al.  (1981). 
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for  different  detection  on  Pinedale  test  tape 
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The  MARS  detector  was  also  one  of  the  better  detectors  although,  on 
the  average,  worse  than  those  mentioned  above.  This  detector  is 
discussed  in  detail  in  Farrell  et  al.  (1981).  The  envelopes  of  several 
narrow  bandpasses  are  computed  and  the  mean  and  standard  deviation  of 
the  relative  maxima  of  these  envelopes  are  exponentially  summed  over  the 
previous  100  maxima.  If  a  fixed  number  of  the  maxima  in  a  2.8  second 
time  window  exceeds  the  mean  plus  a  constant  times  the  standard  devia¬ 
tion,  then  a  detection  is  declared.  This  is  a  "voting"  detector  as 
defined  by  Wirth  et  al.  (1976),  and  for  this  reason,  one  is  not 
surprised  to  find  that  it  does  not  have  the  best  performance. 

The  "online"  detector  is  the  detector  most  like  the  detector 
programmed  for  LASA  by  IBM  and  is  in  routine  use  at  the  SRC  at  this 
writing.  First,  the  data  is  processed  through  a  fixed  filter  and  a  1.5 
second  short  term  average  (with  1.0  second  overlap)  of  the  absolute  data 
values  is  compared  to  an  LTA  exponentially  averaged  over  60  seconds. 
Three  successive  threshold  crossings  give  a  detection.  After  detection, 
the  LTA  update  time  is  switched  to  6  seconds  until  the  signal  shuts  off. 

A  more  detailed  discussion  of  the  actual  imp lemen ta t ion  of  all  of 
the  above  detectors,  together  with  condensed  program  listings,  may  be 
found  in  Appendix  I. 


A  possible  objection  to  the  study  by  Blandford  et  al.  (1981)  is  that 
the  data  were  not  "real",  that  is,  the  phase  of  the  noise  had  been 
scrambled  to  ensure  that  there  were  no  signals  in  the  "noise".  By  so 
doing,  the  authors  may  have  altered  the  statistics  of  the  noise  (for 
example,  by  making  it  more  Gaussian)  in  such  a  way  that  a  detector  which 
was  best  on  that  data  would  not  be  best  in  practice. 

To  test  this  hypothesis,  we  first  ran  the  online  detector  on  BFAK 

beams  for  a  time  span  for  which  experienced  analysts  had  made  especially 

careful  analyses  of  the  develocorder  beam  traces.  As  discussed  in 
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Appendix  I,  the  result  of  this  was  that  the  detector  seemed  to  be 
out-performing  the  analysts,  so  that  there  was  no  way  in  which  the 
analysts  could  function  as  a  standard  against  which  the  different 
detectors  could  be  compared. 

The  remainder  of  this  report  discusses  the  results  of  using  the 
performance  of  the  detectors  on  a  single  instrument  as  a  measure  of 
detector  quality,  while  using  the  analyst  detections  of  small 
teleseismic  signals  on  the  beam  as  a  standard. 
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EXPERIMENTAL  DESIGN 

We  examined  the  A-list  for  November  1,  2,  4,  5,  6  and  found  the  40 
lowest  amplitude  signals  measured  at  BFAK  which  were  confirmed  as  P 
arrivals  from  events  at  distances  of  greater  than  20°.  We  then  searched 
the  bulletins  and  the  raw  analyst  records  and  found  two  4-hour  periods 
during  which  no  arrivals  were  detected. 

To  set  the  false  alarm  rate,  we  first  ran  the  online  detector  on  the 
direct  sum  trace  and  on  a  single  high-gain  instrument  for  the  two 
4-hour  periods.  Any  detection  on  the  single  trace,  which  is  not 
classified  as  a  spike,  local,  or  drop-out,  which  is  classified  as  a 
teleseism,  and  which  is  not  detected  on  the  beam  trace,  is  defined  to  be 
a  false  alarm.  The  threshold  was  adjusted  until  the  false  alarm  rate 
(FAR)  was  in  the  range  3-4/hour.  (Examination  of  Figures  1  and  2  shows 
that  variation  of  the  FAR  over  this  range  corresponds  to  a  variation  in 
magnitude  threshold  of  less  than  0.01  m^ . ) 

We  then  ran  the  other  detectors  on  the  single  channel  data.  If  a 
detection  occurred  within  30  seconds  of  a  spike,  local,  or  drop-out  as 
declared  by  the  online  detector,  then  it  was  not  classified  as  a  false 
alarm,  nor  was  it  a  false  alarm  if  a  detection  occurred  within  30 
seconds  of  a  detection  on  the  beam.  Again,  the  threshold  was  adjusted 
for  each  detector  until  the  false  alarm  rate  was  3-4/hour. 

It  is  of  interest  to  remark  that  during  these  quiet  8  hours  when 
simultaneous  (often  to  the  tenth  of  a  second)  detections  would  occur  on 
the  beam  and  single  instrument,  then  the  beam  S/N  value  would  be  sub¬ 
stantially  higher.  Also,  although  all  the  parameters  were  set  iden¬ 
tically,  the  detection  rate  was  2-4  times  higher  on  the  beam  than  on  the 
single  instrument.  These  facts  certainly  suggest  that  many  real  signals 
were  being  detected  by  the  computer  which  the  analyst  did  not  detect. 

The  40  signals  mentioned  in  the  first  paragraph  of  this  section  were 
analyzed  by  starting  the  detectors  15  minutes  before  the  signal  start 
time  and  seeing  if  a  signal  was  detected  between  the  times  10  seconds 
before  to  30  seconds  after  the  analyst's  start  time.  (The  precursor  10 
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seconds  is  to  allow  for  analyst  error  and  for  possible  window  overlap 
effects  in  the  advanced  detectors.) 

In  the  process  of  recovering  these  40  signals  from  the  station  log 
tapes,  there  were  difficulties  with  6  tapes,  leaving  34  signals  to  be 
analyzed.  Upon  examining  the  detection  results,  we  found  that  the  Z 
detector  often  missed  signals  if  there  was  another  signal  detected  by  it 
a  few  minutes  before,  and  that  the  Walsh  detector  sometimes  failed  if  it 
did  not  have  a  full  13  minutes  of  noise  to  fill  the  detection  histogram 
before  the  signal  arrival. 

The  reason  for  the  Z  detector's  problem  is  that  the  LTA  decay  rate 
was  set  to  be  2  minutes,  and  this  led  to  an  elevated  LTA  after  events, 
and  thus  to  a  missed  signal.  A  similar  problem  may  exist  for  the 
present  coding  of  the  MARS  detector. 

The  Z  problem  can  perhaps  be  solved  by  adjusting  the  LTA  decay  rate 
to  a  shorter  time  for  a  minute  or  so  after  every  detection;  the  Walsh 
detector  weakness  should  not  be  a  problem  in  operational  practice. 

To  eliminate  these  problems,  we  further  discarded  2  events  whose  log 
tapes  began  5  and  7  minutes  before  the  signal,  and  also  eliminated  the 
17  events  with  Z  detections  in  the  5  minutes  before  the  predicted  signal 
when  the  Z  FAR  was  set  to  1/hour.  This  left  a  total  of  only  15  events 
out  of  the  original  40.  While  this  is  a  severe  reduction  in  the  data 
base,  it  seemed  to  be  the  only  way  to  evaluate  the  algorithms  fairly  for 
their  detection  potential  alone.  Any  of  the  detectors  can  and  must,  in 
practice,  be  surrounded  with  substantial  "special  case"  software  to 
handle  the  complicated  situations  which  arise  when  multiple  weak, 
moderate,  and  strong  signals  arrive,  but  the  detector's  success  in  these 
cases  is  little  indication  of  the  underlying  quality  in  pure  detection. 
It  is,  of  course,  in  the  handling  of  these  complex  situations  where  the 
human  analyst's  capabilities  are  most  impressive. 
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RESULTS  AND  DISCUSSION 


In  the  analysis  of  the  15  single  instrument  recordings,  the  Z 
detector  detected  all  15,  the  online  detector  detected  14,  the  Walsh 
detector  detected  13,  and  the  MARS  detector  detected  11.  This  ranking 
is  in  general  agreement  with  the  relative  performance  found  by  Blandford 
et  al.  (1981),  except  that  the  performance  of  the  online  detector  is 
much  better  in  the  present  instance.  This  may  be  due  simply  to  a 
statistical  fluctuation  in  such  a  small  sample,  or  it  may  be  that  taking 
the  absolute  amplitude  is  a  good  thing  to  do  in  real  data  as  suggested 
by  Goforth  and  Herrin  (1981),  Finally,  it  is  most  likely  that  the 
standard  detector  was  at  a  disadvantage  in  the  earlier  study  because  the 
signal  time  window  was  held  at  2.5  seconds  while  the  signal  windows  of 
the  other  detector  were  much  larger  as  appropriate  to  the  larger  signals 
which  were  buried  in  noise. 


Since  this  is  a  new  data  source,  it  might  be  said  that  none  of  the 
detectors  have  been  fairly  evaluated  since  none  of  them  have  been 
optimized  across  parameter  space.  I  have  pretty  much  set  the  parameters 
at  values  which  were  close  to  those  used  previously,  modified  by  what 
would  seem  reasonable  for  the  present  application.  However,  my  feeling 
is  that,  since  these  results  are  roughly  consistent  with  the  performance 
found  by  the  Blandford  et  al.  (1981),  the  referenced  results  can  be 
trusted  as  a  guide  to  practice. 

To  the  extent  that  this  is  true,  then  either  the  online  detector 
modified  to  square  the  data  and  detect  on  one  threshold  crossing,  or  the 
Walsh  detector  would  be  the  "best  buy"  since  they  are  both  very  fast, 
being  limited  by  the  recursive  filter  calculation.  My  own  preference 
would  be  for  the  modified  online  detector  because  it  follows  naturally 
from  classical  detection  theory.  If  there  is  no  problem  with  computer 
time  (the  optimum  Z  detector  runs  a  single  channel  16  times  faster  than 
real  time  on  the  11/70  under  UNIX  as  a  sole  user,  a  dedicated  11/70 
analyzing  several  channels  would  run  2-3  times  faster  per  channel)  then 
the  optimum  detector  would  be  the  obvious  choice.  It  optimally  adapts 
to  changes  in  the  noise  spectrum  including  changes  in  the  microseisms 
and  thus  requires  less  "site  tuning"  and  it  is  the  most  sensitive. 


It  is  extremely  important  to  note  that  before  any  detector  is  relied 
upon  to  select  segments  to  return  to  a  central  site  for  further 
analysis,  it  is  important  that  it  be  evaluated  as  to  its  capability  to 
reject  false  alarms,  to  not  miss  later  phases  and  mixed  events, and  to 
not  "over-produce"  in  the  coda  of  an  event.  To  my  knowledge,  no 
detector  is  currently  functional  which  meets  all  these  criteria, 
although  I  believe  the  online  system  at  the  SRC  is  close,  being  weak 
mostly  in  the  criteria  of  later  phases  and  mixed  events. 

Another  interesting  result  of  this  study  is  that  we  again  see  how 
difficult  it  is  to  find  a  standard  by  which  to  judge  detectors.  Had  the 
detectors  been  only  slightly  better,  or  the  array  only  slightly  smaller, 
then  all  detectors  would  probably  have  detected  all  15  events  and  it 
would  have  been  impossible  to  rank  them.  Possibly  the  only  way  to  get  a 
good  sample  of  real  data  for  comparing  enhancements  of  the  optimum 
detector  will  be  not  from  analyst  picks  but  by  relying  on  the  detector 
operating  on  the  beam  as  a  flag  for  candidate  signals  on  the  single 
instruments.  This,  however,  is  "flying  blind"  and  careful  thought  would 
be  required  to  ensure  that  errors  are  not  being  made. 
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APPENDIX  I 

April  7,  1981  Geotech  Memorandum  from  R.  R.  Blandford, 
"Evaluation  of  SDAC  On-Line  Detector  on  Nov.  1,  2,  1980 
On-Line  ALK  Data  -  Revision  of  30  March  1981  Memo 


TO:  R.  Van  Nostrand 

FROM:  R.  R.  Blandford,  R.  Baumstark 

DATE:  April  7,  1981 

SUBJECT:  Evaluation  of  SDAC  On-Line  Detector  on  Nov.  1,  2,  1980 
On-Line  ALK  Data  -  Revision  of  30  March  1981  memo. 

cc:  M.  Shore,  P.  Kovacs,  A.  Hill,  Capt.  P.  Terry,  Maj  .  G.  Ullrich, 

D.  von  Seggern,  Z .  Der,  R.  Romine 


Summary 

We  have  run  the  on-line  detector  on  the  Nov.  1  and  2,  1980  data  base 
which  was  picked  by  experienced  analysts.  For  speed  of  analysis,  we  used 
the  DP  tapes  as  a  data  source.  In  terms  of  initial  teleseisms,  the  detec¬ 
tor  works  very  well,  detecting  50/57  *  88%  of  all  signals  and  42/44  =  95% 
of  those  not  designated  as  "possible  background",  "probable  background" 
or  "possible  dram.  In  addition,  the  on-line  detector  detected  14  signals 
in  the  first  4  hours  of  Nov.  1  which  were  not  detected  by  analysts  but 
which  upon  close  examination  with  regional  beams  and  filters  we  felt 
were  good  teleseismic  picks.  The  14  signals  are  labeled  A,B,C.  Both 
analysts  felt  that  the  "A"  signals  speak  for  themselves  without  beams  or 
filters;  that  the  B  signals  were  good  and  the  C  signals  questionable. 

One  of  the  two  analysts  felt  that  beams  and  filters  did  not  actually 
held  him  to  decide  if  any  signals  were  real  but  did  offer  more  informa¬ 
tion  about  signals  assumed  to  be  real.  It  is  important  to  state  that  in 
our  opinion  these  are  judgement  calls  and  there  is  no  way  we  know  to  be 
certain  that  one  plot  is  signal  and  the  other  is  noise. 

At  the  14  in  four  hour  rate,  there  would  be  143  missed  by  the  experi¬ 
enced  analyst  in  the  41  hours  of  DPS  tape  data  available.  There  is  no 
particular  reason  to  think  that  these  are  not  "real  signals."  This  is 
not  to  say  that  the  Geotech  analysts  are  "better"  than  experienced  ana¬ 
lysts,  only  that  they  could  concentrate  more  because  they  did  not  have  to 
do  the  routine  scanning. 

In  Table  I  we  see  a  calculation  to  show  that  in  one  day  we  would 
expec t  there  to  be  46  signals  on  a  monitored  channel  (46  x  5  *  230  over 
all  of  ALK)  that  would  be  "unconf irmable"  at  ALK  in  the  sense  that  the 
signal  would  not  be  automatically  detected  at  3  or  more  stations  in 
total.  (The  calculation  assumes  that  there  are  5  m^  =  4.0  events  per 
day  on  the  average,  and  that  the  exponential  law  for  seismicity  holds 
with  b  =  -1.0.  The  S/N  ratio  is  assumed  to  vary  across  the  network  with  o 
*  0.4,  characteristic  of  LASA.)  The  number  of  uncon f i rmab 1 e  detections 
would  be  79  in  41  hours.  Note  that  with  higher  or  lower  false  alarm 
rates,  detection  numbers  could  change  by  a  factor  of  two  and  there  would 
be  no  way  to  "confirm"  that  a  weak  detection  was  a  true  signal  or  false 
alarm.  Note  that  most  of  the  "unconfirmed"  detections  come  from 
magnitude  3.0  events.  This  calculation  is  of  course  only  for 
plausibility.  Thresholds  might  be  lowered  when  looking  for  confirming 
events,  but  in  this  case,  even  more  "one-station"  detections  will  be 
made,  further  complicating  the  issues.  However,  using  the 
machine-Geotech  analyst  as  a  standard  then  the  probability  of  d^.’ction 
of  the  experienced  analyst  is  only  *  50/(50  +  143)  “  25%.  Clearly  this 


makes  it  very  difficult  to  use  the  experienced  analyst  as  a  standard. 
Although  the  experienced  analysis  is  only  "preliminary"  even  a  two-fold 
increase  in  detections  will  leave  the  situation  extremely  ambiguous  and 
impossible  to  interpret  satisfactorily. 

If  there  is  no  standard,  then  it  seems  impossible  to  quantitatively 
compare  one  detector  to  another  using  on-line  data.  (.Except  that  if  a 
detector  was  markedly  inferior  to  the  present  on-line  system  it  could  be 
found  to  be  so.  For  example  a  careful  study  of  all  the  data  by  experienced 
analysts  may  make  it  possible  to  set  thresholds  such  that  the  automatic 
detectors  detect  all  of  the  first  pass  detections  of  the  experienced 
analysts  plus  other  good  picks  without  any  questionable  picks.  This  would 
show  that  the  DP  is  superior  to  the  analyst  in  this  narrow  field.  But  the 
results  of  TR-81-7  would  suggest  that  it  would  be  impossible  to  assert 
that  weak  detections  by  the  best  automatic  detectors  are  or  are  not 
signals  since  most  analyst  "false  alarms"  in  that  paper  looked  like 
signals  and  many  true  signals  detected  by  the  automatic  detector  showed 
little  evidence  of  signal.  The  same  phenomena  are  noted  in  our  on-line 
system  where  the  detector  triggers  on  arrivals  which  would  not  be 
selected  by  the  analyst  unless  he  had  knowledge  from  other  stations 
as  to  the  signal's  shape  and  arrival  time.) 

For  these  reasons,  we  suggest  that  the  aims  (comparing  teleseismic 
detectors)  of  this  detection  study  be  rethought.  If  it  is  truly  desired 
to  compare  different  detectors  then  perhaps  more  "scrambled  phase"  data 
of  different  types  should  be  created.  Another  line  of  work  would  be  to 
run  two  detectors  simultaneously  (higher  frequency,  longer  window  for 
the  second  detector)  to  improve  detections  of  locals  which  were  only 
detected  with  -  50%  reliability,  and  to  test  possible  improvements  in 

detection  of  later  phases  by  adjusting  the  LTA  update  time  constant  and 
threshold  in  the  first  30  seconds  to  1  minute  after  an  initial 
detection. 

Teleseismic  Detection  Results 

The  on-line  detector  at  the  SDAC  is  exactly  as  discussed  in  detail 
in  an  Appendix  to  the  recent  report  by  Blandford,  Racine  and  Romine 
(1981,  SDAC-TR-8 1-7 ) .  Briefly,  a  1.2-3  Hz  prefilter  is  followed  by  an 
STA/LTA  detector  with  an  STA  of  1.5  seconds  and  a  recursively  updated 
LTA  of  ~  30  seconds.  The  STA  is  shifted  0.5  seconds  at  a  time  and  3 
threshold  crossings  in  a  row  are  required  for  detection.  In  running  the 
on-line  detector  on  DPS  tapes  the  time  period  1324  to  2045  Nov.  2  had  to 
be  skipped  because  of  data  gaps. 

A  12  March  1981  VSC  memo  from  R.  Kimmel  transmitted  experienced 
analysis  results  for  01  and  02  Nov,  1980  of  an  infinite  velocity  beam 
from  BFAK.  The  analysis  broke  detections  down  into  event  groups 
with  initial  and  later  phase,  and  the  events  characterized  as  dram 
(locals),  and  pint  or  quart  (multiple  or  impulsive  teleseism).  Remarks 
such  as  possible  background,  possible  dram  (of  a  pint),  possible  noise, 
probable  noise,  possible  regional  motion  (of  a  pint),  possible  dram 
motion  (of  a  quart),  and  possible  surface  motion  (of  a  pint)  indicate 
analyst  uncertainty  over  some  picks.  There  are  a  total  of  13  such  picks 
and  44  with  no  uncertainty  indicated  for  a  total  of  57.  Of  the  57,  50 


were  detected  by  the  machine  for  a  probability  of  detection  of  88%.  If 
we  consider  only  the  44  "more  certain”  picks,  then  42  of  them  were 
detected  for  95%  probability  of  detection.  These  statistics  are 
tabulated  in  Table  II. 

In  Figure  1,  we  see  plots  of  the  two  events  missed  by  the  machine 
and  about  which  the  analyst  indicated  no  uncertainty.  The  upper  signal 

is  rather  high  frequency - perhaps  a  regional  event  although  not  called 

that  by  the  analyst.  Perhaps  also  it  is  simply  the  temporary 
interference  of  the  microseisms  revealing  the  on-going  high  frequency 
noise--that  is--a  false  alarm.  The  lower  event  is  certainly  weak, 
although  we  agree  that  it  should  be  picked.  In  Figure  2  we  see  plots  of 
the  5  signals  not  detected  by  the  machine  out  of  the  population  of  13 
teleseisms  of  which  the  analysts  were  uncertain.  The  second  and  third 
signals  from  the  top  do  show  indications  of  being  "rolling”  background. 

In  Figure  3,  we  see  14  signals  detected  by  machine  in  the  first  4 
hours  of  Nov.  2  which  were  not  picked  by  the  experienced  analyst  but 
which  seem  to  us  to  be  every  bit  as  "real”  as  the  analyst  picks.  In  al¬ 
most  every  case,  there  is  a  definite  change  in  waveform  character  a round 
the  indicated  detection  time  which  persists  well  beyond  the  2.5  second 
window  needed  for  detection  triggering.  As  discussed  in  the  summary, 
this  leads  to  an  estimate  of  analyst  probability  of  detection  of  only 
25%  if  the  machine-analyst  were  adopted  as  a  standard.  This  can  not  be 
thought  to  be  impossible  when  one  considers  the  recent  results  of  the 
analyst-detector  comparison  on  signals  buried  in  phase-scrambled  noise 
in  SDAC-TR-8 1-7  . 

In  Figure  4  we  see  two  signals  which  were  discarded  by  the  analyst 
on  the  basis  that  they  were  not  thought  to  be  signals.  Again,  one 
analyst  used  beams  and  filtering  to  make  this  judgement,  the  other  did 
not  find  these  aids  to  be  a  help  in  the  narrow  question  of  pure 
detection,  although  he  felt  that  they  are  a  great  help  in  putting  out  a 
bulletin.  Beam  plots  for  all  signals  are  available  for  inspection.  It 
seems  worth  remarking  that  Figure  4  may  show  valid  signals.  Good 
detections  were  made  in  TR-81-7  that  looked  as  poor,  and  the  on-line 
detector  has  triggered  on  signals  as  poor  which  have  been  shown  by 
analysis  of  the  complete  event  to  be  real. 


Detection  Timing 

As  discussed  in  the  Appendix  of  SDAC-TR-8 1 -7  ,  the  on-line  detector 
has  a  detection-time  refinement  feature  in  which  the  STA/LTA  is  "backed 
up”  using  a  lower  threshold  than  that  used  to  determine  detection.  It 
has  been  our  experience  that  this  detector  is  generally  better  at 
picking  weak  starts  than  is  the  human  analyst  as  judged  by  its  better 
agreement  with  start  times  determined  by  overlay  techniques  using 
stronger  signals  from  nearby  stations.  This  may  be  seen  to  some  degree 
in  Table  III  where  the  median  machine-analyst  time  is  -0.3  sec  (machin 
late)  for  signals  of  which  the  analyst  was  certain,  but  about  which  the 
analyst  indicated  some  uncertainty  over  start  time  or  signal  type.  The 
early  detection  for  weak  as  opposed  to  strong  signals  suggests  that  the 
automatic  detector  is  more  accurate.  A  definitive  test  would  have  to 


involve  signals  buried  in  noise  or  perhaps  a  formalization  of  the  analyst 
comments. 

Later  Phases 


One  of  the  most  spectacular  failings  of  the  present  detector  is  its 
failure  to  detect  later  phases.  The  present  on-line  system  very 
rapidly  builds  up  the  LTA  after  a  detection  and  then  returns  to  the  very 
slow  LTA  rate  (30  seconds)  existing  before  the  detection.  In  Table  I  we 
see  that  only  29/74  later  phases  were  detected.  The  reasons  for  the 
present  LTA  post-detection  update  parameters  are  strictly  historical  and 
we  are  free  to  change  them.  Perhaps  a  change  for  1  minute  after  an 
initial  detection  to  a  10-sec  LTA  integration  time  coupled  with  a 
lowered  threshold  will  result  in  better  later  phase  detection. 

Local  Signals 

The  on-line  system  has  been  designed  to  detect  small,  impulsive, 
short-period  teleseismic  signals.  The  filter  design  assumes  t*  =  0.3. 
For  regional  phases  0.1  would  be  better  and  a  signal  window  of  5-20 
seconds  would  seem  appropriate.  We  would  run  both  detectors  in  parallel 
and  see  if  a  more  useful  assortment  of  detections  is  made  and  whether  AA 
will  then  turn  out  a  better  Alaskan  regional  bulletin. 
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Figure  3-1.  Machine  detections  missed  by  experienced  analysts  and  rated 
good  picks  by  Geotech  analysts.  First  4  hours  of  Nov.  1. 
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Figure  3-2.  Machine  detections  missed  by  experienced  analysts  and  rated 
good  picks  by  Geotech  analysts.  First  4  hours  of  Nov.  1. 
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Figure  3-3.  Machine  detections  missed  by  experienced  analysts  and  rated 
good  picks  by  Geotech  analysts.  First  4  hours  of  Nov,  1. 


Figure  l* .  Machine  picks  rejected  bv  Oeotech  analysts  as  not 
showing  good  signal . 
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A 


Theoretical  calculation  of  confirmed  and  unconfirmed  detections  at  network  of  5  equal  stations. 
P(D  “  i),  probability  of  exactly  i  out  of  5  detecting.  N(i)  number  of  i-station  events/day. 
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TAPLE  TIT 


Time  Residuals  of  Teleseismic  Analyst  Detections  Also  Detected  by  Machine 
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APPENDIX  II 

Program  Descriptions  and  Condensed  Listings 


Z-De  tec  tor 


The  Z-detector  is  basically  the  same  as  the  standard  detector  except  that 
an  optimum  S/N**2  filter  is  used  instead  of  a  fixed  filter,  the  STA  is 
computed  as  the  log  of  the  sum  of  the  squared  filtered  data  instead  as  of  the 
sum  of  the  absolute  values,  and  finally  detection  is  on  the  Z  variable  equal 

to  (STA  -  LTA)/o.  where  a  is  the  standard  deviation  of  the  STA  about  the  LTA. 

After  this  introduction,  we  see  the  condensed  listing  of  the  routine 
which  replaces  STA/LTA  for  the  standard  detector.  Condensation  is 
indicated  by  three  dots  in  a  vertical  line.  The  omitted  code  can 
generally  be  found  by  reference  to  the  complete  listing  of  the  STA/LTA 
routine  for  the  on-line  detector  which  can  also  be  found  in  this 
Appendix.  A  few  lines  of  code  such  as  type  declarations  unique  to  the  Z 
detector  are  missing  from  the  listing;  for  a  full  understanding  of  the 
code,  the  original  listings  should  be  consulted.  A  full  discussion  of 
the  standard  detector  can  be  found  in  Blandford  et  al.  (1981). 

In  the  calling  sequence,  'cigs'  is  the  standard  deviation  which  is  passed 
back  to  the  calling  routine,  'dpsb',  as  are  the  LTA,  QR ,  etc.  The  variables 

'spec',  'rsp',  and  ' tsa'  are  respectively:  the  noise  spectrum  updated  to  use 

in  the  optimum  filter,  the  amplitude  spectrum  of  the  instrument  response,  and 
the  spectral  shape  due  to  t*.  The  "np"  variables  have  to  do  with  the  length 
of  the  FFT  window  used  to  compute  the  optimum  filter.  In  our  case  np=256,  and 
np2=129,  the  number  of  raw  data  points  to  be  filtered  and  the  dimension  of  the 
spectral  arrays  brought  along. 

First,  the  tapers  are  created,  then  the  initial  values  are  created,  using 
either  the  standard  recursive  filter  or  the  optimum  filter.  We  shall  assume 
throughout  the  following  that  the  optimum  filter  is  used. 

Subroutine  'cfil',  which  may  be  found  following  'STA/LTA',  actually 
applies  the  filter  to  the  data  and  is  quite  straightforward.  The  filter 
is  applied  to  the  first  256  points  of  the  RAW  data.  The  result  is 
placed  in  the  'filter'  array,  the  76th  point  of  which  is  equivalenced  to 
the  first  point  of  the  filter  array.  Since  the  'STA'  calculation  loop 
(Do  1000)  runs  over  NPTS*150  points  beginning  at  filter  (1)  it  runs  only 
to  point  75+150*225  of  the  256  filtered  points  and  so  there  is  no 
problem  about  running  over  the  end  of  the  filtered  data. 


Subroutine  'Rfine'  has  been  modified  to  work  with  the  Z  variable  instead 
of  the  STA/LTA  ratio;  and  thus  it  is  necessary  to  pass  to  it  the  'cig*  value. 

The  'LTA1  and  'cig'  values  are  computed  recursively,  just  as  the  'LTA' 
value  is  updated  in  the  standard  dector,  except  that  we  do  not  allow  the  'cig1 
to  be  updated  when  a  large  signal  arrives,  since  the  transition  from  signal  to 
niose  is  not  characteristic  of  the  variance  of  either  signal  or  noise  by 
itself. 

After  the  Do  1000  loop,  we  have  the  update  of  the  noise  spectrum  using  the 
parameters  set  at  the  beginning  of  the  routine.  The  equation  'tsud=1200'  sets 
the  exponential  updating  of  the  spectrum  to  1200  seconds,  or  20  minutes, 
'dtsud'  is  set  to  15,  in  accordance  with  the  fact  that  the  spectra  are 
updated  only  when  a  new  set  of  15  seconds  worth  of  data  has  been  read  in. 
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Subroutine  Stalta  (CHANID,  RAW,  ICHAK,  CHAN,  XX,STA, 

&  at atim, c igs , LTA,  QR,  DETN, DETECT, STATUS, 

&  spec, rap, taa, irf ,np,np2,np2d) 

Parameter  (ISR  ■  10) 

Parameter  (S  -  20,  NS  »  3,  QP  =  3) 
c...  STALTA  computes  the  short  &  long  term  averages  of  the 

c...  filtered  input*  z“(Lta-sta)/cig  is  compared  to 

c...  a  threshold  value  to  define  a  detection, 
c...  filter  values  squared  and  Log  taken 

c 

c...  Initial  Setup  Section 


do  1  i”l ,np 

1  rawb(i)“1.0 

call  taper(rawb,np, .2) 

npt*np/10 

do  2  i*l,npt 

2  tap( i)“rawb( i) 


c 

c  spectral  update  time;should  eventually  be  put  in  control 

c  parameters 

tsud”1200. 

dtsud“15. 

e2sudKl *0-exp(-dtsud/tsud) 
sigsud^l *0-e2sud 

c 

c...  Compute  initial  LTA  and  thresholds 
c 

if ( irf *eq.l)Then 

Call  Rfil  (NPTS,BCOEF(0,ICHAN(12)),RAW(151). filter) 
else 

Do  17  i“l,np 

17  rawb( i)*RAW( i) 

Call  Detmd(ravb,np,0) 

Call  Taper(ravb,np, .1) 

Do  18  i-l,np 

rawc( i)»Cmplx(rawb( i) ,0 .0) 

18  continue 

Call  Nlogn(8,rawc,-1.0) 
do  19  i“l,np2 

19  spec( i)“2.*(Real(rawc( i) )**2+Aimag(rawc( i) )**2) 

cal  1  cf il( RAW, np, spec .specr ,np2,f ilteq,rsp, tsa,npt, tap) 

end  if 

LTA  -  0 

cig-0. 

cig8“0 . 

Do  20  i  «  1 ,S*NS 
if (irf *ne.l)then 
LTA-LTA+(f ilter( i+90) )**2 
else 

LTA*LTA+(f ilter( i) )**2 
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end  if 

20  continue 

LTA-Log( LTA/ (Float( S*NS) ) ) 

If  (first)  Then 
first  *  .false. 

Cal  1  Header 
End  if 

saqout  E  .false, 
statim  =  timei+dtsta 
do  21  i-l,S*NS 
if(irf  .ne.Dthen 

c  igs  =c  igs  +( LTA-Log ((filter( i+90) )**2) )**2 
else 

cigsEcig8+(LTA-Log( ( f ilter( i) )**2) )**2 
end  if 

21  continue 
cigs=cigs/Float( S*NS) 

DETECT(4)  *  0.0 
DETECT( 5)  -  0.0 

Call  Report  (CHANID, statim, 'Initial  ',0.0, LTA, 
&  statim, saqout, DETECT(4) ,curzr) 

End  if 
c 

c...  Get  channel  calibration  and  delay 
c 


isud*l 

c 

c  isud=l  if  spectrum  to  be  updated 

c 

if(  irf.eq.Dthen 
c 

c...  Load  up  previous  filtered  data  to  avoid  start-up  ring 
c...  Filter  the  data  to  be  used  in  the  STA/LTA  calculation 
c 

Do  40  i  E  1,75 
40  filteq(i)  «  XX(i) 

Call  Rfil  (NPTS,BCOEF(0,ICHAN(12) ),RAW(76) .filter) 

Do  42  i  E  1,75 

42  XX( i)  =  f ilteq( i+150) 

else 

call  cf il( RAW, np, spec. specr .np2, f ilteq.rsp, tsa.npt, tap) 
end  if 
c 

c...  Set  up  display  times  and  files 


nsta  *  1 

ena-0 . 

cig»8qrt(cigs) 

Do  1000  filptr  -  l.NPTS.S 

cursta  “  0.0 

il  -  f ilptr-( S*(NS-1) ) 
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i2  -  f ilptr+(S-l) 

Do  100  i  *  il » i2 

100  cursta  ■  cursta+(f ilter( i))**2 

cursta^LogCcursta/ (Float(S*NS) )) 
if(cig.eq.O . )cig*0 .1 
curz“( cursta-LTA)/cig 
nsta  “  nsta+1 


Do  110  i  -  l.QP-l 
STA(i)  -  STA(  i+1) 

110  QR(i)  -  QR(  i+1) 

STA(QP)  *  cursta 
big  ■  .false. 


Call  Refine  (f ilteq.rawptr,  STA(QP-l),  LTA,  detptr.cig) 


curzr«=log((DETECT(2)-DETECT(3)  )/cig) 

Call  Report  (CHANID .statim, STATUS, 

&  DETECT(2) ,DETECT(3) ,DETECT(1) .saqout,DETECT(4) .curzr) 

isud'O 
c 

c...  if  detect  don't  update  spectrum 

c 

End  if 
c 

c...  Update  LTA  every  'R'  STA's; 
c 

150  Continue 

If  (nsta.eq.R)  Then 
If  (DETN)  Then 
If  (valid)  Then 

LTA  ■  e2sigd*cursta  +  sigdf*LTA 

c  no  cig  update  if  signal  start  where  STA/LTA 

c  very  large 

Else 

valid  *  .true. 

End  if 
Else 

LTA  ■  e2sig  *cursta  +  sigf  *LTA 

cigs-  e2sig  *(LTA-cursta)**2+sigf*cigs 

cig-sqrt(cigs) 

End  if 

tsrt  ■  TRNON 
tend  -  TRNOFF 
nsta  “  0 
ena»ena+l . 

End  if 

c...  Update  start  and  end  detection  thresholds, 
c 


th  ■  tsrt 
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If  (DETN)  th  =  tend 
Continue 

if(isud.eq.l)then 
ncnt=ncnt+l 
do  999  i«i,np2 

spec( i)*e2sud*s pecr( i)+sigsud*spec( i) 

continue 

end  if 

End 
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subroutine  cf il(RAW,np, spec, specr,nP2,f ilteq, rap, tea, npt, tap) 
complex  f iltex(256) ,cf 

"  mJ^RAWC*)  »spec(*)  ,specr(*)  ,filteq(*)  ,rsp(*)  .t»a(*>  .enp.tapC*) 
integer*4  i,np,np2,npt 
do  10  i“l,np 

10  filteq(i)“RAW(i) 

call  detrndCf ilteq, np,0) 
c  call  taper(f ilteq, np, .1) 

do  11  i“l,npt 

filteq(i)=filteq(i)*tap(i) 

11  filteq(np-i+l)“f ilteq(np-i+l)*tap(i) 

do  20  i=l >np 

20  filtex(i)“cmplx(f ilteq(i) ,0.0) 
call  Nlogn(8,filtex,-1.0) 
do  25  i=l  ,np2 

25  specr (i) “Real (filt ex( i) )**2+Aimag( f lit ex( i) )**2 

f iltex(l)“cmplx(0. ,0.) 

do  30  i“2 ,np2 

if ( spec (i) . eq. 0. ) spec (i)“8pec( 128) 

cf=rsp(i)*tsa(i)/8pec(i) 

filtex(i)=filtex(i)*cf 

if (i.ne.nP2)f iltex(np-i+2)*f iltex(np-i+2)*cf 

30  continue 

call  Nlogn(8,f iltex,l .0) 
enp=np 

do  40  i“l ,np 

40  f ilteq (i)“enp*Real(filtex(i) ) 
return 
end 
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Walsh  Detector 

The  basic  idea  of  the  Walsh  detector  is  that,  first,  the  data  are 
frequency  filtered  as  in  the  standard  detector  to  remove  the  micro¬ 
seisms,  then  selected  Walsh  coefficients  are  calculated  for  successive 
32  point  windows  of  10  sps  data.  The  successive  windows  overlap  by  16 
points;  the  Walsh  coefficients  are  prewhitened  with  a  recursively  up¬ 
dated  noise  spectrum,  which  is  updated  only  when  there  are  no  threshold 
crossings.  The  sum  of  the  absolute  values  of  the  whitened  coefficients 
is  calculated  and  compared  to  a  threshold  which,  if  crossed  for  two 
successive  windows,  results  in  a  detection. 

One  of  the  most  interesting  features  of  the  Walsh  detector  is  the 
method  of  setting  the  threshold.  The  successive  'STA'  values  are  used 
to  maintain  a  histogram  and  the  threshold  is  calculated  as  the  sum  of 
the  histogram  median  plus  'xk'  times  the  difference  between  the  75% 
point  and  the  median.  1 xk '  is  typically  in  the  range  3. 0-3. 7.  This 
procedure  for  setting  the  threshold  is  non-parametr ic  and  there  is 
evidence  that  it  is  an  excellent  procedure  which  should  be  considered 
for  any  detection.  it  need  not  be  restricted  in  application  to  the 
Walsh  detector. 

Moving  to  a  consideration  of  the  code,  thee  parameters  S=16  and  NS=2 
provide  for  a  32  point  window  overlapping  by  16  points,  QP=2  provides 
for  the  requirement  of  2  successive  detection  (Q/QP)  where  Q  is  set 
equal  to  2  in  the  "data  cards".  'indlo'  and  'indhi'  are  the  "indif" 
bounding  indices  of  the  Walsh  'sequences'  whose  amplitudes  are  to  be 
calculated  for  detection.  'bakgln'  is  the  number  of  overlapping  32 
point  spectral  windows  used  to  compute  the  detection  histograms.  In  the 
actual  calculations  reported  here  this  value  is  set  to  512,  amounting  to 
an  'averaging'  time  of  13.6  minutes.  ' xk '  is  the  threshold  mentioned 
above . 

The  subroutine  'cwalrt'  computes  the  'indif'  Walsh  functions  ' i y ' 
which  are  of  length  32.  The  routine  'walspc'  computes  the  coefficients 
of  'a'  of  these  Walsh  functions  weighted  by  the  coefficients  ' swt ' . 

After  the  initialization  the  'swt'  values  are  updated  recursively,  near 
the  end  of  the  subroutine. 
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Until  enough  time  has  passed  so  that  the  required  number  of  values 
have  been  computed  to  fill  the  histogram  (nh.gt.nhp),  the  detection  is 
determined  using  the  recursively  updated  LTA  calculated  from  the  STA 
values  going  into  the  histogram.  After  this  time  the  criteria  is 
(curs t a . gt . hid) ,  and  'hid',  the  threshold,  is  calculated  as  discussed 
above . 
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Subroutine  St alt a  (CHANID,  RAW,  ICHAN,  CHAN,  XX,STA, 

&  statim.LTA,  QR,  DETN, DETECT, STATUS) 

Parameter  (ISR  *  10) 

Parameter  (S  “16,  NS  *  2,  QP  *  2 ,indlo»8,indhi*21 , 

&  bakgln*512,indif*14) 
c 

c  Waleh  Transform  version 

c  version  for  single  channel,  for  more  channels  arrays 

c  must  be  defined  in  dpsw  to  save  val, time, hid, swt, and 

c  other  variables,  and  must  be  passed  through  calling 

c  sequence, 

c 


real  swt(indif) ,hld ,a( indif ) .val(bakgln) 
real  time(bakgln) 

integer*2  iy(S*NS, indif)  ,i50,i75,nsta 


c  Initial  set  up 


do  7  i*l .bakgln 

val( i)*i 
7  time(i)*i 
nh-0 

c  setting  histogram  variables  to  some  initial  state 

c 


nhp**bakgln 
i50*(bakgln+l) /2 
i75*( i50+bakgln+l ) /2 
i*0 

do  8  k*indlo , indhi 
j  j-k-1 
i*i+l 
mm*  5 

call  cwalrt(mm, j j,iy(l ,i) ) 

8  continue 

c 

c...  Compute  initial  LTA  and  thresholds 
c 

Call  Rfil  (NPTS , BCOEF(0 , ICHAN ( 1 2) ) ,RAW( 1 51 ) .filter) 

nn-S*NS 

nnn*nn 

nnn2*nn 

do  21  i*l , indif 
21  swt(i)*l .0 

call  walspc(f ilter.nnn, indif ,iy ,a, swt , curs ta) 
do  20  i*l, indif 
20  swt(i)*abs(a(i)) 
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curstaBindif 

LTA«*indif 


Call  Report  (CHANID,statim,'Initial  ',0.0,LTA, 
&  statim, saqout .DETECT (4) ,val(i50) ) 

xk-TRMON 
th“2.2 


Do  40  i  ■  1,75 
40  filteq(i)  B  XX(i) 

Call  Rfil  (NPTS,BCOEF(0,ICHAH(12)),RAW(76) .filter) 
Do  42  i  ■  1 ,75 

42  XX(i)  B  £ilteq( i+150) 


nata  "  1 

Do  1000  filptr  -  1,NPTS,S 
curata  “  0.0 

11  -  £ilptr-(S*(NS-l) ) 

12  =  £ilptr+(S-l) 

if  (  i2  .gt  .NPTS)  il“NPTS-mm2 

avnf«75 

nnn“nn 

call  wal8pc(£ilter(il) ,nnn,indi£,iy, a, swt, curata) 

nata  B  nata+l 

statim  B  statim  +  dtsta 


Do  110  i  -  l.QP-l 

If  ( ( cur a t a. ge . th*LTA. and . ( nh . le . nhp ) ) . or . ( (nh . gt . 

|  &  nhp) .and.cur8ta.gt .hid))  Then 

big  ■  .true. 

i  : 

Call  Report  (CHANID, statim, STATUS, 

&  DETECT(2) .DETECT(3) ,DETECT(1) .aaqout ,DETECT(4) ,i50) 


150  Continue 
c 

c  update  histogram  with  cursta,  and  weights  swt  also 
c 

if(big.eq. .false. ) then 
nhhp*nhp 

call  hi8t(val, time, nhhp, cursta) 
hld«val(i50)+xk*(val(i75)-val(i50) ) 
nhMnh+l 

do  200  iBl,indif 
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200 


8wt(i)*e2sig*abs(a(i) )+sigf*swt(i) 
end  if 

If  (nsta«eq.R)  Then 


End  if 

tsrt  -  TRNON 
tend  -  TRNOFF 
nsta  B  0 


End 
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MARS  Detector 

The  general  flow  in  this  program  is  similar  to  that  in  the  Z-detector 
program  in  that  blocks  of  256  data  points  are  filtered  and  then  analyzed  from 
point  76  to  point  225  filteq.  First  the  envelopes,  'e',  are  calculated  fr  the 
256  point  window,  using  filters  with  center  frequencies  'fc',  and  Q  values 
qe=6*fc.  We  use  9  windows  which  span  the  frequency  range  from  .875  to  3.125 
Hz  using  the  same  spacing,  0.25  Hz,  as  used  at  NORSAR  by  Farrell  et  al  (1981). 

The  Q  values,  qe=6*fc  are  different  from  those  of  12*fc  used  by  Farrell  et 
al  in  order  to  give  the  same  overlap  of  the  bands  in  frequency  as  those 
authors  for  the  successful  performance  of  the  MARS  detector  on  the  PWY  data 
where  the  Q's  were  also  12*fc  but  the  filters  were  spaced  at  only  0.125  Hz. 
qe=6*fc  also  gives  filter  response  times  of  about  6  seconds,  the  probable 
signal  length  for  this  data  set  which  contains  signals  from  smaller  events 
than  those  on  the  VSC  detection  tapes. 

After  the  envelopes  are  computed  using  subroutine  'nbe',  the  envelope 
maxima  in  a  running  window  of  2.7  seconds  are  calculated  and  it  is  determined 
if  they  lie  above  a  threshold.  (This  window  length  is  the  "dispersion 
parameter"  in  the  notation  of  Farrell  et  al;  their  value  for  this  parameter 
was  2.8  seconds.  In  the  present  study  it  is  2.7  secondss  because  of  the 
method  of  creating  this  window.  The  threshold  is  parameterized  as  the 
amplitude  parameter  y  — in  the  program  ,gam'--and  is  set  equal  to  1.8,  the 
same  value  as  in  Farrell  et  al.  The  actual  threshold  for  each  band  is  equal 
to  the  mean  of  the  maxima  for  that  band  plus  gam*  (standard  deviation  about 
the  mean  of  the  maxima).  The  mean  and  standard  deviation  Cnue'  and  'cige'  in 
the  program)  are  calculated  via  recursive  exponentiation  such  that  the 
e-folding  time  is  100x0.9=90  seconds. 

The  "bandwidth"  parameter  of  Farrell  et  al  was  set  equal  to  1  Hz  (k=4)  for 
NORSAR  and  0.62  Hz  (k=5)  for  PWY.  In  this  study  we  select  0.75  Hz  (k=3).  The 
k  parameter  is  signified  in  the  present  program  by  'kd'. 

In  the  program  the  parameter  'nfc2'  is  the  index  of  the  filter  whose 
envelope  values  are  printed  out  in  the  event  of  a  detection  and  has  no  other 
function  in  the  calculations. 
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Subroutine  Stalta  (CHANID,  RAW,  ICHAN,  CHAN,  XX.STA, 

&  statim.LTA,  QR,  DETN .DETECT, STATUS) 

Parameter  (ISR  =  10) 

Parameter  (S  m9,  NS  15  3  ,  QP  ■  3  ,np»256  ,nf  c*9  ,kd“3  ,nfc2*,6) 
c 

c  version  for  MARS,  one  channel  only 

c 


Data  fc/3., 2. 75, 2. 5, 2. 25, 2., 1.75, 1.5, 1.25,1./ 

c 

c...  initial  setup  section 
c 

gam=TRN0N 
do  11  i=l,nfc 
11  qe(i)=  6 .*fc(i) 

df f=float(ISR/2)/float(np/2) 
al=log(2.)/2. 

alp=sqrt(al/3 .14159)  j 

do  10  i=l ,nfc  I 

do  10  j=l,np/2+l  j 

f=f loat( j)*df f 

ex=al*( (qe(i)/fc(i) )**2)*(f-fc(i) )**2 
if (ex.gt .20. )ex=20. 

f il(i, j)=alp*exp(-ex)  ! 

10  continue  " 

nps=np  j 

nfc8=nfc 

call  nbe(RAW,nps,e,nfcs,f il,nps/2+l)  j 

do  9  i=l ,nfc  j 

mue(i)=0. 
do  8  j=2 ,np-l 

if (e(i, j ) .gt .e(i, j-1) .and .e(i, j) .gt .e(i, j+1) 

&  .and.e(i,  j)  .gt .mue(i))mue(i)*e(i,  j) 

8  continue 
mue(i)“mue(i)/l .8 

cige( i)=0 .52*mue( i)  j 

9  continue 
LTA=mue(nf c2) 


Call  Report  (CHANID .statim, 'Initial  '.0.0.LTA, 

&  statim, saqout ,DETECT( 4) ,cige(nfc2) ) 


c  envelope  calculation 

cal  1  nbe(RAW,nps,e,nfcs,f il,nps/2+l) 


npkt*100. 
op  8=1 ./npkt 
lmeps“l .-eps 
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Do  1000  filptr  -  l.NPTS.S 
cursta  “  0.0 

11  «  f ilptr-(S*(NS-l) ) 

12  ■  filptr+(S-l) 
nk-0 

Do  100  i  *  1 ,nfc 
ngam(i)=0 

do  100  j“75+il+l ,75+i2-l 

if(e(i, j) .gt.e(i, j-1) .and.e(i, j).gt. 

&  e(i, j+1) .and .ngam( i) .eq.0)then 
if (i.eq.nfc2)cursta“e(i, j) 
if(e(i, j) .gt.(mue(i)+gam*cige(i) ))then 
nk=nk+l 
ngam(i)“l 
end  if 

mue ( i ) mep s*e ( i , j ) + lmep a*mue ( i ) 
if(i.eq.nfc2)LTA"mue(i) 

cige( i)“8qrt(eps*(mue(i)-e(i, j ) )**2+lmeps*cige(i)**2) 
end  if 

100  continue 


If  (nk.ge.kd)  Then 
nks"nk 
big  =  .true. 


Call  Report  ( CHANID, s tatim, STATUS, 

&  DETECT (2) ,DETECT(3) ,DETECT(1) ,saqout,DETECT(4) ,f loat(nks)) 


End 
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subroutine  nbe(x ,np , e ,ne , f il ,mp ) 
mteger*4  ne ,  i  ,np ,mp ,  ii,  ifd 

real  x(*) ,e(ne,np) ,fil(ne,np/2+l)  ,y(256)  ,s(129)  , 
&z(256) 
do  110  i=l  ,np 
110  y(i)=x(i) 

call  taper(y,np, .2) 
do  100  ii=l,ne 
do  90  i=l ,mp 
90  s(i)=f il(ii,i) 
do  80  i=l ,np 
80  z(i)=y(i) 
ifd=2 

call  fmphl( z ,np, s ,mp , ifd) 
c 

c  have  phaseless  filter  trace 

c 

do  70  i=l ,np 
70  e(ii,i)=z(i) 
ifd=3 

call  fmphl(z,np,s ,mp,ifd) 
c 

c  z  is  now  the  hilbert  transform 

c 

do  60  i=l ,np 

60  e( ii, i)=sqrt(e( iiii)**2+z( i)**2) 

100  continue 
return 
end 
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c 

c...  Detection  Criterion  Common 


c 

c . . . 
c  •  •  • 
c .  .  . 
c .  .  . 
c .  .  . 
c. .  . 
c .  .  . 
c .  .  . 
c .  .  . 
c .  .  . 
c 

c .  .  . 
c  •  «  • 
c .  .  . 
c 

c .  .  . 
c .  .  . 
c .  .  . 
c .  .  . 
c  •  •  • 
c  •  •  • 


NPTS  -  No.  of  DATA  points  to  process 

S  -  No.  of  DATA  points  to  average  for  intermediate  sta  (Y) 

NS  -  No.  of  Y's  to  average  for  STA 

R  -  Interval  in  STA  units  between  LTA  update 

Q  -  q  in  q/q'  criterion 

QP  -  q'  in  q/q'  criterion 

TRNON  -  Turn-on  Threshold 

TRNOFF  -  Turn-off  Threshold 

TSUBC  -  Normal  LTA  Decay  Time  constant 

TSUBCD  -  Decay  Time  constant  during  detection 


IDX  -  Coherency  distance  for  Spkchk  routine 

BRATIO  -  'Big  Spike'  ratio  "  "  " 

SRATIO  -  'Small  Spike'  ratio  "  "  " 


IFLO 

IFMID 

IFHI 

IFCUT 

SPECT 

BCOEF 


-  Low  FFT_channel  for  Spectral  Ratio  Calculation 

-  Middle  “  "  "  "  "  " 

-  High  "  "  "  "  " 

-  Cutoff  FFT_channel  for  inverse  filtering 

-  Local/Teleseism  Discrimination  Ratio 

-  B-coef ficients  for  recursive  filters 


Real  TRNON,  TRNOFF,  TSUBC,  TSUBCD 
Real  BRATIO,  SRATIO,  SPECT,  BCOEF(0:6,3) 
Integer*2  IDX,  NPTS,  R,  Q 
Integer*2  IFLO, IFMID, IFHI ,  IFCUT 
Common  /rit/  TRNON,  TRNOFF,  TSUBC,  TSUBCD, 
&  BRATIO,  SRATIO,  SPECT,  BCOEF, 

&  NPTS,  R,  Q,  IDX,  IFLO, IFMID, IFHI,  IFCUT 
c 

c  End  of  rit.h 
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c...  Time  of  Day  variables 


c 


c  •  •  • 

TIMEI 

-  Start  time 

of 

current  'nsec'  buffer 

c  •  •  • 

TIMERS 

-  Start  time 

of 

run 

c  •  •  • 

TIMERF 

-  Finish  " 

it 

ii 

c .  •  • 

TIMETS 

-  Start  time 

of 

tape 

c  •  •  • 

TIMETF 

-  Finish  " 

n 

II 

Real  TIMEI 
Real  TIMERS,  TIMERF 
Real  TIMETS,  TIMETF 
Integer  IYR,IDOY,IHR,IMN,ISC,ICS 
Integer  SYR.SDOY 
Integer  FYR.FDOY 
c 

c...  I/O  File  Definition  variables 


Integer  TAPFD,  LCNT 
Character  DRIVNO 
Character  DBNAME*7 
Integer*2  RUNSRT(5) 

Integer  SAQLUN ,  EPXLUN 
Character* 15  SAQFIL,  EPXFIL 
Character*80  COMMIT 
Character*?  SMLSTR 

Common  /run/  TIMEI, TIMERS, TIMERF, TIMETS, TIMETF, 
&  IYR, IDOY ,IHR, IMN , ISC , ICS , 

&  SYR.SDOY,  FYR.FDOY, 

&  SAQLUN,  EPXLUN, 

&  LCNT,  RUNSRT, 

&  TAPFD,  DBNAME,  COMMNT ,  SMLSTR, 

&  DRIVNO,  SAQFIL, EPXFIL 
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Subroutine  Stalta  (CHANID,  RAW,  ICHAN,  CHAN,  XX.STA, 

&  LTA,  QR,  DETN, DETECT, STATUS) 

Parameter  (ISR  “  10) 

Parameter  (S  “  5,  NS  =  3,  QP  =  3) 
c 

c...  STALTA  computes  the  short  &  long  term  averages  of  the 

c...  filtered  input.  The  ratio  of  STA/LTA  is  compared  to 

c...  a  threshold  value  to  define  a  detection, 
c 

Include  'rit.h' 

Include  'run.h' 
c 

Character* 14  CHANID 
Integer*4  ICHAN(*) 

Real*4  CHAN(*) 

Real  RAW(*) ,  XX(*) ,  STA(*) ,  LTA 
Integer*2  QR(*) 

Logical  DETN 
Real  DETECT(6) 

Character  STATU S*10 
c 

Real  calib,  delay 
Complex  xcl(64),  xc2(64) 

Real  dt,  df ,  dtlta,  dtsta,  fltlta 
Real  qsum,  sumlo,  sumhi,  frnax 
Real  e28ig,  sigf,  e2sigd,  sigdf 
Real  statim,  rawtim.f iltim 
Real  xrl(64),  xr2(64) 

Real  filteq(300) ,  filter(150) 

Equivalence  (f ilter(l) ,f ilteq(76) ) 

Real  tsrt, tend ,th,  cursta 
Integer*2  i,j,  il,i2 
Integer*2  iknt,  nsta 

Integer*2  detptr ,f ilptr.rawptr , f ftptr 
Logical  big,  qqp,  first 
Logical  detsrt,  defend,  saqout 
Logical  drop,  BSPK,  SSPK, valid 
c 

Data  first/. true./ 

Data  filteq  /300*0.0/ 
c 

c...  Initial  Setup  Section 

c...  Done  only  if  LTA  *  -999.9 

c 

If  (LTA. eq. -999. 9)  Then 
c 

c...  Define  scaling  factors 

c 

dt  -  1 .0/Float( ISR) 

df  -  1 .0/ (64*dt ) 

dtsta  “  S*dt 

dtlta  -  (R*S)*dt 

e2sig  ■  1.0  -  Exp(-dtlta/tsubc) 

sigf  “1.0-  e2sig 

e2sigd  “  1.0  -  Exp(-dtlta/tsubcd) 

sigdf  “1.0-  e2sigd 
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Apr  8 

c 

c  •  •  • 
c 


20 


& 

c 

c  •  •  « 

c 

c 

c .  .  . 
c .  .  . 
c 

40 

42 

c 

c  •  •  • 
c  •  •  • 
c 


c 

c .  .  . 
c 


c 

c  •  •  • 
c  •  •  • 
c 
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Compute  initial  LTA  and  thresholds 

Call  Rfil  (NPTS,BCOEF(0,ICHAN(12) ),RAW(151) , filter) 
LTA  =  0 

Do  20  i  =  1 ,  S*NS 

LTA  -  LTA+Abs (f ilt er( i) ) 

If  (first)  Then 
first  =  .false. 

Call  Header 
End  if 

saqout  -  .false, 
statim  =  timei+dtsta 
fltlta  -  Float(LTA) /(S*NS) 

DETBKT(4)  =0.0 
DETECT(5)  =0.0 

Call  Report  (CHANID , statim, 'Initial  ' ,0 .0 , fltlta, 
statim, saqout ,DETECT(4) ) 

End  if 

Get  channel  calibration  and  delay 

calib  =  chan(5) 
delay  =  chan(6) 

Load  up  previous  filtered  data  to  avoid  start-up  ring 
Filter  the  data  to  be  used  in  the  STA/LTA  calculation 

Do  40  i  =  1,75 

fiiteq(i)  =  XX(i) 

Call  Rfil  (NPTS,BC0EF(0,ICHAN(12)),RAW(76), filter) 

Do  42  i  =  1,75 

XX(i)  =  f ilteq(i+150) 

Set  up  display  times  and  files 
Filter  time  delay  is  0.3  sec  for  Rfil 

rawtim  =  timei-15.0+delay 
filtim  =  timei-15 .0+delay-0 .3 
statim  ■  timei-  7  .5+delay-0 .3-dtsta 

Define  thresholds 

tsrt  -  TRN0N*LTA 
tend  =  TRNOFF*LTA 
th  =  tsrt 

If  (DETN)  th  =  tend 

Set  up  STA  &  LTA  Storage 
Begin  Loop  over  data 

nsta  “  1 

Do  1000  filptr  »  1  ,NPTS , S 
cursta  ■  0.0 

11  =  f ilptr-(S*(NS-l) ) 

12  ■  filptr+(S-l) 
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Do  100  i  *  il,i2 

cursta  “  cursta+Abs(filter(i)) 
nsta  *  nsta+l 
statim  *  statim  +  dtsta 

Shift  detection  flag  array; 

Check  for  STA/LTA  >  threshold; 

Do  110  i  «  1,QP-1 
STACi)  ■  STA( i+1) 

QR(i)  ■=  QR(i+l) 

STA(QP)  =  cursta 
big  *  .false. 

QR(QP)  *  0 

If  (cursta. ge.th)  Then 
big  ■  .true. 

QR(QP)  -  1 
End  if 
qsum  *  0 
Do  112  i  *  1 , QP 

qsum  “  qsuai+QR(i) 
qqp  ■  .falser 

If  (qsum.ge.Q)  qqp  *  .true. 

Check  for  detection: 

big  =  true  if  last  STA  >  threshold 

qqp  **  true  if  q/ q'  STA's  >  threshold 

DETN  *  true  if  a  detection  is  in  progress 


detsrt 

detend 


.false. 

.false. 


Detection  Start? 

If  (big  .and.  .not. DETN  .and.  qqp)  Then 
detsrt  *  .true, 
valid  ■  .true. 

DETN  *  .true, 
th  ■  tend 

rawptr  *  filptr+75 

Call  Refine  (filteq, rawptr,  STA(QP-l),  LTA,  detptr) 
DETECT(l)  -  filtim  +  (detptr-l)*dt 
DETECT (2)  -  0.0 
Do  120  i  -  1,QP 

If  (STA(i).gt.DETECT(2))  DETECT(2)  -  STA(i) 
Continue 

DETECT(2)  -  DETECT(2)/Float(S*NS) 

DETECT(3)  -  LTA/Float(R*S) 

Check  for  detection  on  data  drop-out. 

iknt  *  0 
drop  ■  .false. 

Do  125  i  -  (detptr-(ISR/4)),detptr+ISR 
If  (RAW(i).ne.O)  Then 


I 


t 


Else 

iknt  ■  iknt+1 
End  if 

If  ( iknt.ge.ISR/2)  drop  *  .true. 

125  Continue 

If  (drop)  Then 

STATUS  =  'Drop-out' 

DETECT(4)  -  -999.99 
valid  *  .false. 

Goto  150 

End  if 
c 

c...  Load  up  FFT  buffer 
c 

fftptr  =  detptr-(64/2) 

Do  130  i  •=  1,64 

j  ■  fftptr+(i-l) 
xrl(i)  =  RAW(j) 

130  Continue 

c 

c...  Check  for  Big  spikes 

c 

Call  Spkchk  (64,  xrl ,  IDX,  BRATIO,  BSPK) 

If  (BSPK)  Then 

STATUS  =  'Spike_l ' 
ichan(lO)  *  ichan(10)+l 
valid  *  .false. 

End  if 
c 

c...  Detrend  &  taper  the  data 

c 

Call  Detrnd  (xrl, 64,0) 

Call  Taper  (xrl, 64, 0.1) 
c 

c...  Load  up  complex  array  &  do  FFT 
c 

Do  132  i  «  1,64 

132  xcl(i)  ■  Cmplx(xrKi)  ,0.0) 

Call  Nlogn  (6,xcl,-1.0) 
c 

c...  Load  up  inverse  transform  buffer 

c...  Remove  low  frequency  components  from  FFT  data; 
c...  Do  inverse  transform 

c 

Do  134  i  ■  1 ,  64 
134  xc2(i)  =  xcl(i) 

c 

xc2(l)  -  (0.0, 0.0) 

Do  140  i  “  2.IFCUT 

xc2(i)  -  (0.0, 0.0) 
xc2(64-( i-2) )  -  (0.0, 0.0) 

140  Continue 

Call  Nlogn  (6,xc2,+1.0) 

Pull  out  inverse  xformed  data 


c 

c . . . 
c 
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I 


I 


ft 


142 


144 

c 

c  •  •  • 

c 


146 


148 

c 

c .  •  . 
c  •  •  • 
c  •  •  • 
c 


c 

c  •  •  • 
c 


Do  142  i  *  1 ,64 

xr2(i)  »  Real(xc2( i) ) 

Check  for  Small  spike 

If  (.not.BSPK)  Then 

Call  Spkchk  (64,xr2 ,1024,SRATI0»  SSPK) 

If  (SSPK)  Then 

STATUS  *  'Spike_2' 
ichan(ll)  “  ichan(ll)+l 
valid  =  .false. 

End  if 

End  if 

Compute  the  amplitude  spectra. 

DETECT(5)  -  0.0 
fmax  =0.0 
Do  144  i  =  1,33 

xrl(i)  -  Sqrt(Real(xcl( i) )**2  +  Atmag(xcl( i) )**2) 

If  (xrl( i) .ge . fmax)  Then 
fmax  *  xrl(i) 

DETECT(5)  *  Float(i-l)*df 
End  if 
Continue 

Compute  the  spectral  ratio 

sumlo  =0.0 

Do  146  i  -  IFLO.IFHID 

sumlo  =  sumlo+xrl(i) 
sumhi  *  0.0 

Do  148  i  =  IFMID+1 ,IFHI 

sumhi  *  sumhi+xrl(i) 

If  valid  detection,  use  Spectral  Ratio  for  Local/Teleseismic 
discrimination.  ALL  detections  with  logSR  <  -0.15  are  called 
teleseisms . 

DETECT(4)  *  -999.99 

If  (sumlo. gt. 0.0  .and.  sumhi. gt .0 .0) 

DETECT(4)  *  Logl0( sumhi/ sumlo) 

If  (DETECT(4) .It. -0.15)  valid  -  .true. 

If  (valid)  Then 

If  (DETECT(4) .gt.SPECT)  Then 
STATUS  ■  'Local' 
ichan(8)  ■  ichan(8)+l 
Else 

STATUS  *  'Teleseiem' 
ichan(9)  •  ichan(9)+l 
End  if 

End  if 
End  if 

Continue  Detection? 
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c 

c  •  •  • 
c 


& 

c 

c  •  •  • 
c 

150 


c 

c  •  •  • 
c 

1000 

c 


If  ( DETN  .and.  qqp)  th  ■  tend 

Detection  End? 

If  (.not. big  .and.  DETN  .and.  .not. qqp)  Then 
detend  =  .true. 

DETN  =  .false, 
th  =  tsrt 
saqout  =  .false. 

If  (STATUS. eq. 'Local'  .or.  STATUS. eq. 'Teleseism  )  Then 
saqout  =  .true. 

End  if 

Call  Report  (CHANID.statim, STATUS, 

DETECT(2) ,DETECT(3) ,DETECT(1) , saqout ,DETECT( A) ) 

End  if 

Update  LTA  every  'R'  STA's; 

Continue 

If  (nsta.eq.R)  Then 
If  (DETN)  Then 
If  (valid)  Then 

LTA  =  e2sigd*cursta  +  sigdf*LTA 
Else 

valid  =  .true. 

End  if 
Else 

LTA  =  e2sig  *cursta  +  sigf  *LTA 
End  if 

tsrt  ■  TRN0N*LTA 
tend  -  TRNOFF*LTA 
nsta  =  0 
End  if 

Update  start  and  end  detection  thresholds, 
th  =  tsrt 

If  (DETN)  th  »  te.id 
Continue 

Return 

End 
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